/*
 *  static.cpp
 *  EPPI-0.0
 *
 *  Created by Juan Carlos Garcia Orden on 15/11/07.
 *
 */

#include "static.h"

Static::Static()
:	Analysis(),
	linSolver(),
	nlSolver(),
	numericTangFlag(),
	vTime(2,0),
	tol(),
	maxIter(0)
{}

// ***

Static::Static(ifstream &in, string & na, string &ty, int di, string &gf, string &gp, double &gs,Simulation *pSim)
:	Analysis(na,ty,di,gf,gp,gs,pSim),
	linSolver(),
	nlSolver(),
	maxIter(30),
	tol(1e-10),
	numericTangFlag(),
	vTime(2,0)
{
	string keyword;
	string loadName;					/** Auxiliary string */
	string loadType;					/** Auxiliary string */
	string bodyName;					/** Auxiliary string */
	string auxs;
	string line;
	char* dispFile;
	string dispDofs;
	vector<string> vdispPointNames;
	string dispPointsetName;
	string sysName;
	Process *procPointer;
	int i;
	double dispVal;
	stringstream stnum;
	string snum;
	vector < string > newtonData;
    double value;
	
	while (in >> keyword && keyword.compare("ENDANALYSIS") )
	{
		if(!keyword.compare("LINSOLVER"))
			in >> linSolver;
		else if(!keyword.compare("NLSOLVER"))
		{       
			getline(in,line);
            istringstream iss(line);
            while (iss >> keyword )
            {
				newtonData.push_back(keyword);
			}
                 
            nlSolver = newtonData[0]; //  As usual the user sets nlSolver, with maxIter and tol being set by default.

            if(newtonData.size() == 2)  // The user only sets maxIter or tol and tol or maxIter is set by default
            {
				stnum << newtonData[1];
				stnum >> value;
				if (value>=1)
				{
					maxIter = value;
				}
				else tol = value;
			}
            if(newtonData.size() == 3) // The user sets maxIter and tol
            {
				stnum << newtonData[1];
				stnum >> value;
				if (value>=1)
				{
					maxIter = value;
					stnum.clear();
					stnum << newtonData[2];
					stnum >> value;
					tol = value;									
				}
				else 
				{
					tol = value;
					stnum.clear();
					stnum << newtonData[2];
					stnum >> value;
					maxIter = value;						
				}
			}
		}
		// edit starts here
		else if(!keyword.compare("TIMESPAN")) 
		{
			createTimeVec(in);  
		}
		else if(!keyword.compare("CONSISTENT"))
		{
			integratorType = 3;	// 3 means consistent integration for internal variables.
            cout << "iType= " << integratorType << endl;
		}
		// edit ends here
		else if(!keyword.compare("REMOVELOAD"))
		{
			changeLoads = true;
			/** Read load to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("REMOVESYSTEM"))
		{
			haveToRebuild = true;
			/** Read system to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("REMOVEBODY"))
		{
			haveToRebuild = true;
			/** Read body to be removed.		***** TO BE DONE */
		}
		else if(!keyword.compare("ADDSYSTEM"))
		{
			haveToRebuild = true;
			/** Read system to be added.		***** TO BE DONE */
		}
		else if(!keyword.compare("ADDBODY"))
		{
			haveToRebuild = true;
			/** Read body to be added.			***** TO BE DONE */
		}
		else if(!keyword.compare("GRAVITY"))
		{
			double grav;
			int gravDir;
			in >> grav;
			in >> gravDir;
			// This load has to be sent to all bodies in all processes
			parentSim->setGravity(grav,gravDir);
		}
		else if(!keyword.compare("AUGMENTED"))
		{
			augLagFlag = true;
		}
		else if(!keyword.compare("BEGINLOAD"))
		{
//					cout << "keyword***** is " << keyword << "\n";
			changeLoads = true;		/** Loads change, because a new load is defined next */
			in >> keyword;
			if(!keyword.compare("NAME"))
			{
				in >> loadName;
			} else cout << "*ERROR: keyword NAME is missing or out of place from load \n";
			in >> keyword;
			if(!keyword.compare("TYPE"))
			{
				in >> loadType;
			} else cout << "*ERROR: keyword TYPE is missing os out of place from load \n";
			in >> keyword;
			if(!keyword.compare("BODY"))
			{
				in >> bodyName;
			}
			else cout << "*ERROR: keyword BODY is missing or out of place from load \n";
			sysName = getStringHead(bodyName);
			procPointer = parentSim->getProcFromSys(sysName);	/** Pointer to the process where the body is */
			procPointer->sendLoadToBody(loadType,sysName,bodyName,in,vTime);	/** Send the load to the body */		
		}
		else if(!keyword.compare("BEGINDISP"))
		{
			//vSysNames.clear();
			in >> keyword;
			if(!keyword.compare("NAME"))
			{
				in >> auxs;
			} else cout << "*ERROR: keyword NAME is missing or out of place from load \n";
			in >> keyword;
			if(!keyword.compare("POINTS"))
			{
				getline(in,line);
				//cout << "line is " << line << "\n";
				istringstream iss(line);
				while (iss >> keyword )
				{
					//in >> auxs;
					vdispPointNames.push_back(keyword);
				}
				// cout << "Ha pedido puntos " << endl;
			}
			else if(!keyword.compare("POINTSET"))
			{
				in >> dispPointsetName;
				// cout << "Pide pointset " << endl;
			}
			in >> keyword;
			if(!keyword.compare("DOFS"))
			{
				// cout << "Pide dofs " << endl;
				in >> dispDofs;
			}
			in >> keyword;
			if(!keyword.compare("FILE"))
			{
				//in >> dispFile;
				in >> auxs;
				dispFile=new char[auxs.size()+1];
				dispFile[auxs.size()]=0;
				memcpy(dispFile,auxs.c_str(),auxs.size());				
				// cout << "Ha pedido file " << endl;
				parentSim->createHistory(dispFile);
			}
			else if(!keyword.compare("VALUE"))
			{
				in >> dispVal;
				parentSim->createHistory(dispVal);
				auxs = "value";
				stnum.clear();
				stnum << dispVal;
				stnum >> snum;
				auxs+="-";
				auxs+=snum;
				dispFile=new char[auxs.size()+1];
				dispFile[auxs.size()]=0;
				memcpy(dispFile,auxs.c_str(),auxs.size());
			}

			// Sending pointset and points
			if(vdispPointNames.size()==0)	// Pointset
			{
				sysName = getStringHead(dispPointsetName);
				procPointer = parentSim->getProcFromSys(sysName);
				procPointer->addImposedDispPointSet(dispPointsetName,dispDofs,dispFile);
			}
			else	// Points
			{
				// Loop on points
				for(i=0; i<vdispPointNames.size(); ++i)
				{
					sysName = getStringHead(vdispPointNames[i]);
					procPointer = parentSim->getProcFromSys(sysName);
					procPointer->addImposedDispPoint(vdispPointNames[i],dispDofs,dispFile);
				}
			}
			vdispPointNames.clear();
		}
	}
	if(!nlSolver.compare("newton-numeric")) numericTangFlag=1; else numericTangFlag=0;
//	if(!linSolver.compare("gauss")) lmx::setLinSolverType(0);
//	else if(!linSolver.compare("cg")) lmx::setLinSolverType(2);
}

// ***
/*
Static::Static(string & na, string & ty, int di, string ls, string nls, bool rebuild, bool reset, Process *p)
:	Analysis(na,ty,di,rebuild,reset),
	linSolver(ls),
	nlSolver(nls)
{}
*/
// ***

double Static::calculate()
{
	if (parentSim->lProcess==0) cout << "Starting calculation of analysis " << name << ", static ... \n";
	double time = 0; // edit: t added ... maybe better, if it is a local variable of analysis or static...
	int numStep = 0; // edit: 
	int i;
	int j;
	int numIter = 0;
// 	int maxIter = 30;
// 	double h;
	bool itConverged;
// 	double tol = 1.0e-10;
	double minTol = tol*tol;	// Tolerance to account for linear problems, that needs just one iteration
	double resmod ;	// resmod is a relative residual energy.  force residual * q increment between iterations / initial energy
	double resmodMax;
	double resmodIni;
	bool haveInterResConverged = true;
	int numInterfaceIter;
//	numInterfaceIter = (parentSim->nProc)*(parentSim->nProc)+1;	// Number of interface iterations
//	numInterfaceIter = (parentSim->nProc)+15;
	numInterfaceIter = (parentSim->nProc)-1;
	vector<double> deltaq(parentSim->lProcess->numDofs);
        
	cout << "maxiter = " << maxIter << endl;
	
	parentSim->lProcess->writeToFile("--------------------------");
	parentSim->lProcess->writeToFile("PERFORMING CALCULATION ...");
	parentSim->lProcess->writeToFile("");
	parentSim->lProcess->writeToFile("Number of interface iterations: ", numInterfaceIter);
	parentSim->lProcess->writeToFile("");
	parentSim->lProcess->initializeRes(0,augLagFlag);	// 0 is modeCalc, not used for statics
	parentSim->lProcess->initializeTang(0,false, false, false, implicitFlag);		// false is for velocity projection and lumped mass, not used for statics
        parentSim->lProcess->createConvergedVectors();
	
	if(printIterFlag==1 && parentSim->lProcess->num==0)	parentSim->lProcess->createIterOutFile();	//Just for one process, all have the same number of iterations		
	
	//h = vTime[1]-vTime[0]; // cambiar!!! paso de tiempo variable, constante
        time = vTime[0];

        
	if(printStressFlag==1 || plotStressFlag==1)	parentSim->lProcess->fillMSmooth(time);	// Smoothing matrix, just once calculation at the begining

	timeA = MPI::Wtime();	// *****TIME INFO
	if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("   Preprocess :", timeA-(parentSim->lProcess->tini));

	//*****Sergio
	int plotStep = 0;
	int nShot;

	nShot = (int) (graphicsStep/(vTime[1]-vTime[0]));

	if (nShot < 1 or graphicsStep == 0.0)
	{
			nShot = 1;
	}
	// start edit for quasistatic calculations
	for (numStep = 0; numStep < vTime.size()-1; numStep++) 			// loop over all timesteps
	{
		timeB = MPI::Wtime();	// *****TIME INFO
        
        parentSim->lProcess->setQConv();	// Store converged value of q
		if (numStep==0)	parentSim->lProcess->setQIni();
		
		parentSim->lProcess->update(integratorType);	// Update data at components

		if(printDofsFlag==1) parentSim->lProcess->writeFileDofs(time);			// write dof data for current time
		if(printStressFlag==1 || plotStressFlag==1)
		{
			parentSim->lProcess->calcSmoothStress(time);
			if (printStressFlag==1)	parentSim->lProcess->writeFileStress(time);
		}

		if(printEnerFlag==1)
		{
			parentSim->lProcess->calcEner();
			parentSim->collectEner();
			parentSim->writeFileEner(time);
		}
		
		if (numStep == plotStep)
		{
			parentSim->lProcess->createGraphOutFile(graphicsFormat,graphicsPath,numStep); 	// create graphic output file for current timestep
			parentSim->lProcess->writeFileGraph(graphicsFormat,graphicsPath,numStep,plotDispFlag,plotProcFlag,plotStressFlag);			// write graphic output file for current  timestep

			plotStep = numStep + nShot;
		}

//		parentSim->lProcess->createGraphOutFile(graphicsFormat, graphicsPath, numStep);	// create graphic output file for current timestep
	//	parentSim->lProcess->writeFileGraph(graphicsFormat, graphicsPath ,numStep,plotDispFlag,plotProcFlag,plotStressFlag);			// write graphic output file for current  timestep
		/*
		if(printIterFlag==1)
		{
		}
		*/

		time = vTime[numStep+1];						// get current time step

		if (parentSim->lProcess->num==0) cout << "-------- \n";
		if (parentSim->lProcess->num==0) cout << "--- Timestep " <<  numStep+1   << " of " << vTime.size()-1;	// print current time step
		if (parentSim->lProcess->num==0) cout << "  (" << time << "s of " << vTime.back() << "s)" <<  "\n";	// print current time step in seconds
		
		parentSim->lProcess->writeToFile("   ");
		parentSim->lProcess->writeToFile("----	Timestep ", numStep+1 );		// print current time step to file
		
// edit ends here
		augIter = 0;
		timeC = MPI::Wtime();	// *****TIME INFO
		if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("         Writting post : ", timeC-timeB);
		
		// Lagrange multiplier loop
		do{
			timeC1 = MPI::Wtime();	// *****TIME INFO
			if (augLagFlag==true)
			{
				if (parentSim->lProcess->num==0) cout << "# - Iteration LAMBDA " << augIter << "\n";
				parentSim->lProcess->updateLagMul();	// Update Lagrange multipliers at components
			}
			//timeC2 = MPI::Wtime();	// *****TIME INFO
			//parentSim->lProcess->writeToFileTime("				Update Lagrange multipliers : ", timeC2-timeC1);
			// Newton-Raphson loop
			for(numIter=0; numIter<maxIter; numIter++)
			{
				timeD = MPI::Wtime();	// *****TIME INFO
				
				itConverged = true;
				if (parentSim->lProcess->num==0) cout << "-------- \n";
				if (parentSim->lProcess->num==0) cout << "----NR Iteration " << numIter + 1 << "\n";
				parentSim->lProcess->writeToFile("----NR Iteration ", numIter+1);
                                
                                //To Calculate the pseudo-velocity in static analisys
                                for (i=0; i < ((parentSim->lProcess->q).size()) ; ++i) 
                                {
                                        (parentSim->lProcess->dq)[i] = 1.0 / h * ( (parentSim->lProcess->q)[i] - (parentSim->lProcess->qConv)[i] );
                                }
//                              cout << "dq[" << 3 <<"]= " << parentSim->lProcess->dq[3] << endl;
			
				// Compute residual of isolated processes
				// First 0 means computation for relative residual energy evaluation
				// Second 0.0 means time step (no applicable here)
				// IntegratorType is 0 for statics or 3 for consistent viscoelastic static.
				
// 				if ( time <= 10.0) h = vTime[1] - vTime[0];
// 				else if(time > 10.0 && time < 50.0) h = 1.0;
// 				else if(time >= 50.0) h = 5.0;
				// h = vTime[1] - vTime[0];

				parentSim->lProcess->formIsolatedResidue(time,numericTangFlag,0,integratorType,h);
				// Compute tangent of isolated processes	
				parentSim->lProcess->formIsolatedTangent(time,numericTangFlag);	// 0 means computation for relative residual energy evaluation
				// Reset internal interface contributions
				parentSim->lProcess->resetInterInt();		
				// Prepare tangent of internal dofs

				if(numIter==0 && numStep==0)	parentSim->lProcess->factorizeInternalMatrix(0);
				else parentSim->lProcess->factorizeInternalMatrix(1);

				timeE = MPI::Wtime();	// *****TIME INFO
				if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("                  Form isolated and factorize Kint : ", timeE-timeD);
				
				for(j=0; j<numInterfaceIter;j++)
				// Loop on the interfaces on each process
				{		
					cout << "		Iteration " << j << " at interfaces in process " << parentSim->lProcess->num << "...\n";
			
					// Set external interface contributions
					// true means that KIB, KBB, etc have to be buit from process before adding new values
					if (j==0) parentSim->lProcess->setInterExt(true);
					else parentSim->lProcess->setInterExt(false);
					
					// Set internal interface contributions
					parentSim->lProcess->setInterInt();
				} 
				timeF = MPI::Wtime();	// *****TIME INFO
				if (printTimeFlag==1)	parentSim->lProcess->writeToFileTime("                     Interfaces iteration : ", timeF-timeE);
				
				// Form final residue and tangent of each process
				parentSim->lProcess->formResidue(0);	// 0 means computation for relative residual energy evaluation

				parentSim->lProcess->formTangent();

				timeF1 = MPI::Wtime();	// *****TIME INFO
				if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("                  Form residual and tangent : ", timeF1-timeF);
				
				// Check convergence
				if (parentSim->lProcess->num==0) cout << "	Checking convergence ...\n";
				parentSim->lProcess->writeToFile("");
			
				resmod = 0.0;
				for (i=0; i<deltaq.size(); ++i)
				{
					resmod+= ((parentSim->lProcess)->res)[i] * deltaq[i];
				}
				resmod = fabs(resmod);

				if (numIter==0)	resmod = 2*tol;										// To force the first iteration
				else
				{
					if (resmodIni!=0.0 && resmod>=minTol) resmod/= resmodIni;	// resmod is a relative residual energy
					//				if (resmodIni!=0.0 && resmod>=minTol) resmod = deltaq.norm2();
				}
		
				MPI::COMM_WORLD.Allreduce(&resmod,&resmodMax,1,MPI::DOUBLE,MPI::MAX);
		
				if (resmod > tol)
				{
					cout << "		Process " << rankk << " has NOT converged, Residue = " << resmod << "\n";
					parentSim->lProcess->writeToFile("		Process ",rankk," has NOT converged, Residue = ",resmod);
					itConverged=false;
				}
				else
				{
					cout << "		Process " << rankk << " has converged, Residue = " << resmod << "\n";
					parentSim->lProcess->writeToFile("		Process ",rankk," has converged, Residue = ",resmod);
					if (resmodMax > tol)		itConverged=false;
				}

				timeG = MPI::Wtime();	// *****TIME INFO
				if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("                  Convergence check : ", timeG-timeF1);

				if(itConverged==true)
				{
					timeG = MPI::Wtime();	// *****TIME INFO
					if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("               Newton iteration : ", timeG-timeD);
					break;
				}
				else
				{
					cout << "	Solving and updating process " << parentSim->lProcess->num << "...\n";
					// Solve linear system and update at each process
					if (numIter==0)	parentSim->lProcess->solveUpdate(deltaq,0,time);
					else parentSim->lProcess->solveUpdate(deltaq,1,time);
			
					if(numIter==0)
					{
						resmodIni = 0.0;
						for (i=0; i<deltaq.size(); ++i)
						{
							resmodIni+= ((parentSim->lProcess)->res)[i] * deltaq[i];
						}
						resmodIni = fabs(resmodIni);
					}
					timeG1 = MPI::Wtime();	// *****TIME INFO
					if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("                     Solve and update : ", timeG1-timeG);
					if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("               Newton iteration : ", timeG1-timeD);
				}
				//break;
			}
			// End of Newton loop
			
			timeG2 = MPI::Wtime();	// *****TIME INFO
			
			if(numIter<maxIter && itConverged==true)	// Convergence is achieved
			{		
				parentSim->lProcess->writeToFile(" ------------------ \n");
				parentSim->lProcess->writeToFile("** All processes converged at iteration ",numIter+1);
				parentSim->lProcess->writeToFile(" ------------------");
				parentSim->lProcess->writeToFile("Solution process ",rankk);
				parentSim->lProcess->writeToFile(parentSim->lProcess->q);
			}
			else
			{
				cout << "*ERROR: Convergence in one or more processes not achieved at maxIter = " << maxIter << "\n";
			}
			if(printIterFlag==1 && parentSim->lProcess->num==0)	// Just for one process, all have the same iterations
			{
				parentSim->lProcess->writeFileIter(time,numIter+1);
			}
			if (augLagFlag==true)
			{
				if (augIter==0)	augConvergenceTot=false;
				else 
				{
					augConvergence = true;
					parentSim->lProcess->checkConvLagMul(augConvergence);	// Check convergence of Lagrange multipliers
					MPI::COMM_WORLD.Allreduce(&augConvergence,&augConvergenceTot,1,MPI::BOOL,MPI::LAND);
				}
			}
			augIter+=1;
			timeH = MPI::Wtime();	// *****TIME INFO
			if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("            Augmented Lagrange iteration : ", timeH-timeC1);
		}	while(augLagFlag==true && augConvergenceTot==false && augIter<maxIterLambda);
		// End of Lagrange multiplier loop
		if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("         Augmented Lagrange loop : ", timeH-timeC);
		if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("      Step : ", timeH-timeB);
	}
	// End of time loop

	if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("   Total steps : ", timeH-timeA);

	if(printDofsFlag==1) 
	{
		parentSim->lProcess->writeFileDofs(time);			// write dof data for current time
		parentSim->lProcess->closeFileDofs();			// close file
	}
	if(printEnerFlag==1)
	{
		parentSim->lProcess->calcEner();
		parentSim->collectEner();
		parentSim->writeFileEner(time);
		parentSim->closeFileEner();
	}
	/*
	if(printIterFlag==1 && parentSim->lProcess->num==0)	// Just for one process, all have the same iterations
	{
		parentSim->lProcess->writeFileIter(time,numIter+1);
	}
	*/
	if(printStressFlag==1 || plotStressFlag==1)
	{
		parentSim->lProcess->calcSmoothStress(time);
		if (printStressFlag==1)
		{
			parentSim->lProcess->writeFileStress(time);
			parentSim->lProcess->closeFileStress();
		}
	}
	parentSim->lProcess->createGraphOutFile(graphicsFormat, graphicsPath, numStep);	// create graphic output file for last timestep
	parentSim->lProcess->writeFileGraph(graphicsFormat, graphicsPath,numStep,plotDispFlag,plotProcFlag,plotStressFlag);			// write graphic output file for last  timestep
	parentSim->lProcess->closeGraphOutFile(graphicsFormat,numStep);	// close output graphic file with specific commands (needed by pvd format)
	
	
	parentSim->lProcess->writeToFile(" ------------------ \n");
	parentSim->lProcess->writeToFile("Calculation finished \n");
	
	timeI = MPI::Wtime();	// *****TIME INFO
	if (printTimeFlag==1) parentSim->lProcess->writeToFileTime("   Writting final post : ", timeI-timeH);
	return 0.0;
}

// ***

// edit starts here
void Static::createTimeVec(ifstream &in)
{
		double t_start; 
		double t_end;    
		//double stepsize; 

		vTime.clear(); 
		in >> t_start;		// start time
		//in >> stepsize;		// time increment
		in >> h;		// time increment
		in >> t_end;		// end time

		// cout << t_start << "\t" << stepsize << "\t" << t_end << "\n"; // borrar

/* 
Two possibilities to create the time vector for a given TIMESPAN e.g t = [0.5  1  3]
	1) creates a time vector that does always contain "t_end": vTime = [0.5  1.5  2.5  3]
	2) creates a time vector that uses a constant increment and therefore may not contain "t_end"
		for example here t will create: vTime = [0.5 1.5 2.5 3.5] 
		note that the last element is rounded: e.g [0.5 1 2.9] will create [0.5 1.5 2.5]

	** change the comments to switch between modes	
*/	
	// 1)
		/*
		for (int i = 0; i <= abs(((t_end - t_start)/stepsize)+0.5); i ++) 
		{
			vTime.push_back(t_start+stepsize*i);				// cout << vTime[i] << "\t"; // borrar
		} 
	
		if ((vTime.back()-t_end) > 0) {
			vTime[vTime.size()-1]=t_end;
		}
		else 
		{
			vTime.push_back(t_end);
		}
*/

	// 2)
		/*
		for (int i = 0; i <= fabs(((t_end - t_start)/stepsize)+0.5); i ++) 
		{
			vTime.push_back(t_start+stepsize*i);
//			cout << "vTime[" << i << "]=" << vTime[i] << "\n";
		} 
		*/
		for (int i = 0; i <= fabs(((t_end - t_start)/h)+0.5); i ++) 
		{
			vTime.push_back(t_start+h*i);
		}
}
// edit ends here

// ***

void Static::calcVel ( void )
{
        
}

